In [2]:
using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations, DifferentialEquations
gr();

M3M6: Methods of Mathematical Physics

$$ \def\dashint{{\int\!\!\!\!\!\!-\,}} \def\infdashint{\dashint_{\!\!\!-\infty}^{\,\infty}} \def\D{\,{\rm d}} \def\E{{\rm e}} \def\dx{\D x} \def\dt{\D t} \def\dz{\D z} \def\C{{\mathbb C}} \def\R{{\mathbb R}} \def\CC{{\cal C}} \def\HH{{\cal H}} \def\I{{\rm i}} \def\qqqquad{\qquad\qquad} \def\qqand{\qquad\hbox{for}\qquad} \def\qqfor{\qquad\hbox{for}\qquad} \def\qqwhere{\qquad\hbox{where}\qquad} \def\Res_#1{\underset{#1}{\rm Res}}\, \def\sech{{\rm sech}\,} \def\acos{\,{\rm acos}\,} \def\vc#1{{\mathbf #1}} \def\ip<#1,#2>{\left\langle#1,#2\right\rangle} \def\norm#1{\left\|#1\right\|} \def\half{{1 \over 2}} $$

Dr. Sheehan Olver
s.olver@imperial.ac.uk


Website: https://github.com/dlfivefifty/M3M6LectureNotes

Lecture 18: Orthogonal polynomials and singular integrals

This lecture we do the following:

  1. Cauchy transforms of weighted orthogonal polynomials
    • Three-term recurrence and calculation
    • Hilbert transform of weighted orthogonal polynomials
    • Hilbert transform of weighted Chebyshev polynomials

Cauchy transforms of orthogonal polynomials

Given a family of orthogonal polynomials $p_k(x)$ with respect to the weight $w(x)$ on $(a,b)$, we always know it satisfies a three-term recurrence: \begin{align*} x p_0(x) &= a_0 p_0(x) + b_0 p_1(x) \\ x p_k(x) &= c_k p_{k-1}(x) + a_k p_k(x) + b_k p_{k+1}(x) \end{align*}

Consider now the Cauchy transform of the weighted orthogonal polynomial: $$ Q_k(z) := {\cal C}_{[a,b]}[p_k w](z) = {1 \over 2 \pi \I} \int_a^b {p_k(x) w(x) \over x -z} \dx $$

Theorem (Three-term recurrence Cauchy transform of weighted OPs) $C_k(z)$ satisfies the same recurrence relationship as $p_k(x)$ for $k=1,2,\ldots$: \begin{align*} z C_0(z) &= a_0 C_0(z) + b_0 C_1(z) - {1 \over 2 \pi \I} \int_a^b w(x) \dx \\ z C_k(z) &= c_k C_{k-1}(z) + a_k C_k(z) + b_k C_{k+1}(z) \end{align*}

Proof \begin{align*} z C_k(z) &= {1 \over 2 \pi \I} \int_a^b {z p_k(x) w(x) \over x -z} \dx = {1 \over 2 \pi \I} \int_a^b {(z -x) p_k(x) w(x) \over x -z} \dx + \int_a^b {x p_k(x) w(x) \over x -z} \dx \\ &= -{1 \over 2 \pi \I} \int_a^b p_k(x) w(x) \dx + \int_a^b {(c_k p_{k-1}(x) + a_k p_k(x) + b_k p_{k+1}(x) w(x) \over x -z} \dx \\ &= -{1 \over 2 \pi \I} \int_a^b p_k(x) w(x) \dx + c_k C_{k-1}(z) + a_k C_k(z) + b_k C_{k+1}(z) \end{align*} when $k > 0$, the integral term disappears.
⬛️

This gives a convenient way to calculate the Cauchy transforms: if we know $C_0(z) ={\cal C}w(z)$ and $\int_a^b w(x) \dx$, solve the lower triangular system: $$ \begin{pmatrix} 1 \\ a_0-z & b_0 \\ c_1 & a_1-z & b_1 \\ &c_2 & a_2-z & b_2 \\ && c_3 & a_3-z &\ddots\\ &&&\ddots & \ddots \end{pmatrix}\begin{pmatrix}C_0(z) \\C_1(z) \\C_2(z) \\C_3(z) \\ \vdots \end{pmatrix} = \begin{pmatrix}C_0(z) \\{1 \over 2 \pi \I} \int_a^b w(x) \dx \\0 \\0 \\ \vdots \end{pmatrix} $$

Example (Chebyshev Cauchy transform)

Consider the Chebyshev case $w(x) = {1 \over \sqrt{1-x^2}}$, which satisfies $\int_{-1}^1 w(x) \dx = {\pi}$. Recall that $$ C_0(z) ={\cal C}w(z) = { \I \over 2\sqrt{z-1}\sqrt{z+1}} $$ Further, we have \begin{align*} x T_0(x) = T_1(x) \\ x T_k(x) = {T_{k-1}(x) \over 2} + { T_{k+1}(x) \over 2} \end{align*} hence \begin{align*} z C_0(z) = C_1(z) - {1 \over 2 \I} \\ z C_k(z) = {C_{k-1}(z) \over 2} +{C_{k+1}(z) \over 2} . \end{align*} In other words, we want to solve $$ \begin{pmatrix} 1 \\ -z & 1 \\ 1/2 & -z & 1/2 \\ &1/2 & -z & 1/2 \\ && 1/2 & -z &\ddots\\ &&&\ddots & \ddots \end{pmatrix}\begin{pmatrix}C_0(z) \\C_1(z) \\C_2(z) \\C_3(z) \\ \vdots \end{pmatrix} = \begin{pmatrix} { \I \over 2\sqrt{z-1}\sqrt{z+1}} \\{1 \over 2 \I} \\0 \\0 \\ \vdots \end{pmatrix} $$ with forward substitution.


In [5]:
x = Fun()
w = 1/sqrt(1-x^2)
z = 0.1+0.1im

n = 10
L = zeros(ComplexF64,n,n)
L[1,1] = 1
L[2,1] = -z
L[2,2] = 1
for k=3:n
    L[k,k-1] = -z
    L[k,k-2] = L[k,k] = 1/2
end



C = L \ [ im/(2sqrt(z-1)sqrt(z+1)); 1/(2im); zeros(n-2)]


Out[5]:
10-element Array{Complex{Float64},1}:
   0.4999250218677838 + 0.004998750393615982im
 0.049492627147416784 - 0.44950762277386im    
 -0.40012497188352847 - 0.08500174951890464im 
 -0.11251727162034156 + 0.35248227849337344im 
   0.3071250618607855 + 0.13299475089351104im 
  0.14734333381379644 - 0.2644583159425141im  
 -0.22476473190952334 - 0.15641774731925456im 
 -0.16101273073185018 + 0.18822182009675853im 
   0.1549178217438016 + 0.16185956519223624im 
  0.15962438204216325 - 0.12486634270955096im 

In [6]:
T₅ = Fun(Chebyshev(), [zeros(5);1])
cauchy(T₅*w,z) - C[6]


Out[6]:
-5.551115123125783e-17 + 5.551115123125783e-17im

Warning This fails for large $n$ or large $z$:


In [11]:
x = Fun()
w = 1/sqrt(1-x^2)
z = 5+6im

n = 100
L = zeros(ComplexF64,n,n)
L[1,1] = 1
L[2,1] = -z
L[2,2] = 1
for k=3:n
    L[k,k-1] = -z
    L[k,k-2] = L[k,k] = 1/2
end

C = L \ [ im/(2sqrt(z-1)sqrt(z+1)); 1/(2im); zeros(n-2)]

T₂₀ = Fun(Chebyshev(), [zeros(20);1])

C[21], cauchy(T₂₀*w, z)


Out[11]:
(-2.2332625421658917e6 + 521568.0612707699im, 0.0 + 8.834874115176436e-18im)

Get around it by dropping the first row:


In [12]:
L[2:end,1:end-1]


Out[12]:
99×99 Array{Complex{Float64},2}:
 -5.0-6.0im   1.0+0.0im   0.0+0.0im  …   0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.5+0.0im  -5.0-6.0im   0.5+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.5+0.0im  -5.0-6.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.5+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im  …   0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im  …   0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
     ⋮                               ⋱                                    
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im  …   0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im  …   0.5+0.0im   0.0+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im     -5.0-6.0im   0.5+0.0im   0.0+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.5+0.0im  -5.0-6.0im   0.5+0.0im
  0.0+0.0im   0.0+0.0im   0.0+0.0im      0.0+0.0im   0.5+0.0im  -5.0-6.0im

In [13]:
C = L[2:end,1:end-1]\ [1/(2im); zeros(n-2)]

C[6]- cauchy(T₅*w, z)


Out[13]:
-3.072062106493749e-17 + 4.410020129853371e-17im

Hilbert transform of weighted orthogonal polynomials

Now consider the Hilbert transform of weighted orthogonal polynomials: $$ H_k(x) = \HH_{[a,b]}[p_k w](x) = {1 \over \pi} \int_a^b {p_k(t) w(t) \over t-x} \dt $$

Just like Cauchy transforms, the Hilbert transforms have

Corollary (Hilbert transform recurrence) \begin{align*} x H_0(x) &= a_0 H_0(x) + b_0 H_1(x) - {1 \over \pi} \int_a^b w(x) \dx\\ x H_k(x) &= c_k H_{k-1}(x) + a_k H_k(x) + b_k H_{k+1}(x) \end{align*}

Proof Recall $$ \CC^+ f(x) + \CC^- f(x) = -\I \HH f(x) $$ Therefore, we have $$ C_k^+(x) + C_k^-(x) = -\I \HH[w p_k](x) $$ hence we have \begin{align*} x H_0(x) &= \I x (C_0^+(x) + C_0^-(x)) = \I \left[a_0 (C_0^+(x) + C_0^-(x)) + b_0 (C_1^+(x) + C_1^-(x)) -{1 \over \pi \I} \int_a^b w(x) \dx \right]\\ &= a_0 H_0(x) + b_0 H_1(x) - {1 \over \pi} \int_a^b w(x) \dx \end{align*} Other $k$ follows by a similar argument.

⬛️

Example: weighted Chebyshev

For $$ H_k(x) = \int_{-1}^1 {T_k(t) \over (t-x)\sqrt{1-t^2}} \dt $$ The recurrence gives us \begin{align*} x H_0(x) &= H_1(x) -1 \\ x H_k(x) &= {H_{k-1}(x) \over 2} + {H_k(x) \over 2} \\ \end{align*} In this case, we have $H_0(x) = \HH[w](x) = 0 $. Therefore, we can rewrite this recurrence as \begin{align*} H_1(x)&= 1 \\ x H_1(x) &= {H_2(x) \over 2} \\ x H_k(x) &= {H_{k-1}(x) \over 2} + {H_k(x) \over 2} \\ \end{align*} This is precisely the three-term recurrence satisfied by $U_{k-1}$! We therefore have $$ H_k(x) = U_{k-1}(x) $$


In [14]:
x = 0.1

T = Fun(Chebyshev(),[zeros(n);1])
hilbert(w*T,x) - Fun(Ultraspherical(1), [zeros(n-1);1])(x)


Out[14]:
0.0

This gives a very easy way to compute Hilbert transforms: if $$ f(x) = \sum_{k=0}^\infty f_k T_k(x) $$ then $$ \HH\left[{f \over \sqrt{1-\diamond^2}}\right](x) = \sum_{k=0}^\infty f_{k+1} U_k(x) $$